
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### MDPREF #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#





####
## The following function, written by
## William G. Jacoby, performs a
## nonmetric multidimensional preference
## scaling analysis. 
####

###     Requires opscale function 
###     from optiscale package
###
###     Assumes larger data values
###     correspond to more positive 
###     responses
###
###     Within-row variance must be nonzero
###     for all rows of the data matrix
###
###    Arguments:
###
###         xdata  Input data frame or matrix
###
###      standard  If TRUE, data will be standardized
###                within rows. If FALSE, data will
###                be centered, but not standardized
###
###      meas.lev  Measurement level of input data, 
###                1 = nominal (not recommended),
###                2 = ordinal
###
###     meas.proc  Measurement process of input data,
###                1 = discrete
###                2 = continuous
###
###        dimens  Dimensionality of scaling solution
###
###
###   mdpref function returns a list with six elements:
###
###         col.coords  Data frame containing coordinates for
###                     column objects
###
###         row.coords  Data frame containing coordinates for
###                     row objects
###
###   predicted.values  Matrix of model-based predicted values
###
###          os.values  Matrix of optimally-scaled data values
###
###            history  Iteration history
###
###          fitvector  Vector of model fit in increasing 
###                     dimensionality from one to rank of data matrix
###
#####
#####


mdpref <- function (xdata, standard = FALSE,
                    meas.lev = 2, 
                    meas.proc = 1,
                    dimens = 2) {
  
  
  n.bad <- length(apply(xdata, 1, var)[apply(xdata, 1, var) == 0])
  if (n.bad > 0) {
    where.bad <- names(apply(xdata, 1, var)[apply(xdata, 1, var) == 0])  
    print(where.bad)
    stop ("Preceding rows of data matrix have no variance" )
  }
  
  var.names <- colnames(xdata)  
  obs.names <- rownames(xdata)
  
  ###
  ###   Standardize or center (the default)
  ###   the entries within rows of the data 
  ###   and assign as initial matrix of
  ###   optimally-scaled values
  ###
  
  if (standard == TRUE) {
    xdata <- t(apply(xdata, 1, scale))
  }
  
  else {
    center <- function (x) {scale(x, center = TRUE,
                                  scale = FALSE)}
    
    xdata <- t(apply(xdata, 1, center))
  }
  
  
  ###
  ###   Initialize the matrix of optimally-scaled
  ###   values, the previous fit, the iteration
  ###   number, and the improvement in fit over
  ###   the previous iteration
  ###
  
  xdata.os <- xdata
  
  prev.fit.model <- 0
  
  niter <- 0
  
  improve = 1
  
  history <- c()
  
  ###
  ###   Start iterations
  ###
  
  while (improve > .01 & niter <= 25) {
    niter <- niter + 1
    
    ###
    ###   Perform SVD on OS version of data
    ###
    
    decomp <- svd(xdata.os)
    
    ###
    ###   Calculate goodness-of-fit
    ###   for specified dimensionality
    ###   and improvement in fit over previous iteration
    ###
    
    d.sqd <- decomp$d^2
    
    fit.model <- sum(d.sqd[1:dimens]) / sum(d.sqd)
    
    fitvector <- cumsum(d.sqd) / rep(sum(d.sqd), length(d.sqd))
    
    improve <- fit.model - prev.fit.model
    
    ###
    ###   Create iteration history
    ###
    
    if (niter == 1) {
      history <- c(niter, fit.model, improve)
    }
    if (niter > 1)  {
      history <- rbind(history, c(niter, fit.model, improve))
    }
    
    ###
    ###   Obtain terminal points of vectors for 
    ###   row objects, calculate predicted
    ###   ranks, norm the row object vectors to unit length
    ###
    
    row.c1 <- decomp$u[, 1:dimens] %*% diag(decomp$d[1:dimens])
    
    pred.xdata <- row.c1 %*% t(decomp$v[,1:dimens])
    
    sumsqd <- apply(row.c1^2, 1, sum)
    
    root.sumsqd <- (matrix(sumsqd, nrow = length(sumsqd), ncol = 1)) ^ .5
    
    row.coord <- row.c1 / (root.sumsqd %*% matrix(1, nrow = 1, ncol = dimens))
    
    ###
    ###   Obtain new optimally scaled data values,
    ###   using predicted ranks from fitted model
    ###
    
    for (i in 1:nrow(xdata.os)) {
      opped <- opscale(x.qual = xdata[i,],
                       x.quant = pred.xdata[i,],
                       level = meas.lev,
                       process = meas.proc,
                       rescale = T)
      xdata.os[i,] <- opped$os
    }
    
    ###
    ###   Set current fit to previous fit for next iteration
    ###
    
    prev.fit.model <- fit.model
    
  }
  
  ###
  ###   The preceding right bracket 
  ###   ends the program loop
  ###
  
  ###
  ###   Create iteration history matrix
  ###
  
  history <- matrix(data = history, nrow = niter, 
                    ncol = 3)
  
  colnames(history) <- c("Iteration", "Fit", "Improvement")
  
  ###
  ###   Create variable points
  ###
  
  variable.coords <- as.data.frame(decomp$v[, 1:dimens])
  
  rownames(variable.coords) <- var.names
  
  ###
  ###   Create data frame with terminal points
  ###   of row vectors
  ###
  
  row.coord <- as.data.frame(row.coord)
  
  rownames(row.coord) <- obs.names
  
  mdprefout <- list(
    variable.coords, row.coord, pred.xdata, 
    xdata.os, history, fitvector)
  
  names(mdprefout) <- c("col.coords", "row.coords",
                        "predicted.values",
                        "os.values",
                        "history",
                        "fitvector")
  
  mdprefout
  
}
